
#  Replication script for
#  Graham and Huber, "Designing Survey Experiments,"
#  forthcoming in the Handbook of Experimental Methodology

library(meta)

d = read.table("C:/R/metadata.tab", sep = "\t", header = TRUE)

#  meta analysis on post design (standard experiment / no pre measure of outcome)
meta_post = metagen(d$post, d$postse)

#  meta analysis of difference between the post and pre-post designs (i.e., difference between standard and within designs)
meta_diff = metagen(d$prediff, d$prediffse)

#  get SE from meta-analysis of pre-post difference
the_se = meta_diff$seTE.random

#  specify range of potential true differences between the two designs
#  (low end = 0.1 = about 10 percent attenuation relative to baseline of 0.104)
attenuation = seq(from = 0, to = .104, by = .002)

#  for each possible true value, simulate power: 
#  % of times estimate is large enough to reject the null of no effect
#  (i.e., % of times estimate is greater than 1.96 SEs)
power = rep(NA, length(attenuation))
set.seed(0)
for(i in 1:length(attenuation)){
  power[i] = mean(rnorm(10000, attenuation[i], the_se) > 1.96*the_se)   #  do this step analytically?
}

#  summarize results
results = data.frame(
  Effect = meta_post$TE.random,
  Attenuation = attenuation,
  `Percent` = attenuation / meta_post$TE.random,   #  calculate % atenuation: divide by baseline from meta analysis on the post design
  Power_sim = power,
  Power_analytic = 1-pnorm(1.96*the_se, attenuation, sd = the_se)
)

results

######################
#  Figure 2         ##
######################

library(tidyverse)

theme_allplots = function(textsize = 10.5){
  theme_bw() +
    theme(
      text = element_text(family = "Times", size = textsize, color = "gray5"),
      strip.text = element_text(size = textsize),
      axis.title.x = element_text(size = textsize, margin = margin(8, 0, 0, 0)),
      axis.title.y = element_text(size = textsize, margin = margin(0, 8, 0, 0)),
      axis.text.x = element_text(size = textsize-1),
      axis.ticks = element_line(linewidth = .25, color = "gray60"),
      plot.margin = margin(t=0, r=1),
      plot.title = element_text(size = textsize+1.5),
      legend.position = "bottom",
      legend.title = element_blank(),
      panel.spacing.x = unit(.5, "cm"),
      strip.background = element_blank(),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank()
    )
}


g = results |>
  ggplot(
    aes(x = Percent, y = Power_analytic)
  ) +
  geom_hline(aes(yintercept = 0), color = "gray70", lty = 2, linewidth = .25) +
  geom_hline(aes(yintercept = 1), color = "gray70", lty = 2, linewidth = .25) +
  geom_line() +
  theme_allplots() +
  scale_x_continuous(labels = scales::percent) +
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = "Percent Attenuation due to Consistency Bias",
    y = "Power"
  )

ggsave("C:/Users/tuq69844/Dropbox/Graham Huber survey experiments chapter/figure2.pdf", g, "pdf", width = 4, height = 3.25)
